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We demonstrate that substantial progress can be achieved in the study of the phase structure of 4-dimensional 
compact QED by a joint use of hybrid Monte Carlo and multicanonical algorithms, through an efficient parallel 
implementation. This is borne out by the observation of considerable speedup of tunnelling between the metastable 
states, close to the phase transition, on the Wilson line. Our approach leads to a general parallelization scheme 
for the efficient stochastic sampling of systems where (a part of) the Hamiltonian involves the total action or 
energy in each update step. 



1. Introduction 

It is embarrassing that lattice simulations of 
compact QED still have not succeeded to clar- 
ify the order of the phase transition near /3 = 1, 
the existence of which was established in the clas- 
sical paper of Guth §. This is mainly due to 
the failure of standard updating algorithms, like 
metropolis or heatbath to move the system at 
sufficient rate between the observed metastable 
states near its phase transition. Due to super- 
critical slowing down (SCSD) the tunneling rates 
decrease exponentially in L 3 and exclude the use 
of lattices large enough to make contact with the 
thermodynamic limit by finite size scaling tech- 
niques (FSS) §. 

Torrie and Valleau Q| used arbitrary "non- 
physical" sampling distributions in their method, 
termed "umbrella sampling" , to improve the effi- 
ciency of stochastic sampling for situations when 
dynamically nearly disconnected parts of phase 
space occur by biassing the system to frequent 
the dynamically depleted, connecting regions of 
configuration space. Berg and Neuhaus Q ap- 
plied the idea of "umbrella sampling" under the 
name "multi-canonical algorithm" (MUCA) to 
the simulation of a variety of systems exhibiting 
first-order phase transitions. In this procedure, 
the biassing weight of a configuration with action 
S is dynamically adjusted such as to achieve a 



near-constant overall frequency distribution over 
a wide range of S within a single simulation. 

2. Multicanonical Sampling (MUCA) 

"Canonical" Monte Carlo generates a sample of 
field configurations, {4>} according to the Boltz- 
mann weig ht P can (0) - e-P 8 ^. The canonical 
action density which in general exhibits a double 
peak structure at a first-order phase transition, 
can be rewritten as 



—ps 



(1) 



with the spectral density p(S) being independent 
of the inverse temperature p. 
The multicanonical approach aims at generating 
a flat action density 

N UVCA (S,/3 C ) = const., for S' maxl < S < S n ^ 2 , (2) 

in a range of S that covers the double peaks lo- 
cated at S maxl and 5 , max2 - Therefore one modifies 
the sampling weights introducing a multicanoni- 
cal potential V MUCA (S), 



P ( <1\ „-/3cS-Vmuca(S) 

r MUCA V 1 C ! 

r io g iv ca „(s maxl ,/3 c ) 

V MVCA {S) = I logiV ca „(S,/3 c ) 

{ logiV ca „(S max2 ,/3 c ) 



(3) 



S<S max i 
S max i<S<S„ 

s>s max2 



which is constant outside the relevant action 
range. Since V r MUCA (S') is unknown at the begin 
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of the simulation, it is instrumental for MUCA 
to bootstrap from good guesstimates. We shall 
do so by starting from an observed histogram 
of the canonical action density, N caa (S, /3 C ), see 
eq. (0), at the supposed location of the phase 
transition, /3 C . From the action density, we com- 
pute V MU ca(S)- The sampling then proceeds with 
the full MUCA weight, 



-$ a S-V MUCA (S) 
-($ c +$(S))s~a(S) 



(4) 



In order to compute expectation values of observ- 
ables, one has to reweight the resulting action 
density to reconstitute the proper canonical den- 
sity P can (5,/3 C ) ~ Pmuca(S) e^ucA(s). The com _ 
putation of the multicanonical weights requires 
the knowledge of the global and not just the lo- 
cal change in action for each single update step. 
As a consequence, MUCA is not parallelizable for 
local update algorithms. For remedy, we propose 
to utilize the HMC updating procedure. 

3. Hybrid Monte Carlo (HMC) 

In addition to the gauge fields <pfi{x) one intro- 
duces a set of statistically independent canonical 
momenta ir^x), chosen at random according to a 

Gaussian distribution exp(— \-). The action S[<f>] 
is extended to a guidance Hamiltonian 



fJ,.X 



(5) 



Starting with a configuration (0, tt) at time t = Q, 
the system moves through phase space according 
to the equations of motion 

1 dH . dH d 

OTTn 



<9^ 



leading to a proposal configuration ((/>', n') at time 
t = t. Finally this proposal is accepted in a global 
Metropolis step with probability 

Pace = min (l, e~ AH ) , AH=H[<f,',n']-ni4>,n]. (7) 

The equations of motion are integrated numer- 
ically with finite step size At along the trajectory 
from t = up to t — N md At — t. Using the 



leap-frog scheme H as symplectic integrator the 
discretized version of eq. @ fulfills the detailed 
balance condition being time reversible and mea- 
sure preserving. Each integration step approxi- 
mates the correct H with an error of 0(Ai 3 ). 

4. Merging MUCA and HMC 
for Compact QED (MHMC) 

We consider a multicanonical HMC for pure 
4-dimensional U(l) gauge theory with standard 
Wilson action defined as 



S[4>] 



E 
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(M*)) 



0^(X) = 4>^(X) + 4>y{x + fl) - + H) - 4>y(x). 

Eq. |^ suggests to consider an effective action S 
including the "multicanonical potential" Vmuca, 

S = $ c S + V MUCA (Sj c ) = c +P(S))S+a(S).{8) 

We now define MHMC making use of the 
multicanonical potential as a driving term 
within molecular dynamics via the multicanon- 
ical Hamiltonian, TL = |^i 2 +S, inducing an 
additional drift term such that the resulting force 
is given by 



AU {x—0)—sin6 llv (x) 



(9) 



The MHMC is governed by the dynamics un- 
derlying the very two peak structure: U MU ca is 
repelling the system out of the hot (cold) phase 
towards the cold (hot) phase, thus increasing its 
mobility and enhancing flip-flop activity Q. 

5. Tunneling Behaviour 

In order to quantify the efficiency of the 
MHMC, we introduce the average flip time r flip 
defined as the inverse number of the sum of flips 
between the two phases multiplied by the total 
number of trajectories. For reference, we addi- 
tionally measured r flip from the Metropolis algo- 
rithm with reflection steps (MRS) which is con- 
sidered as a very effective local update algorithm 
for U(l) @- 

1 Note that MHMC requires the computation of the 
global action (to adjust the correct multicanonical weight, 
eq. (W)) at each integration step along the trajectory of 
molecular dynamics to guarantee reversibility. 
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^flip 


MHMC 

r fli P 
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1.001600 


508(12) 


650(20) 
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1.007370 


1023(60) 


1173(50) 
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1.009300 


2474(117) 


2006(84) 


12 


1.010143 


5470(770) 


3260(440) 


14 


1.010668 


16400(3300) 


5090(630) 


16 


1.010753 


44800(9700) 


6350(860) 



for the various lattices of volume V = L 4 . 



With the results for r aip on lattices up to 16 4 
we are in the position to estimate the scaling 
behaviour of MHMC in comparison to standard 
MRS updates. According to the expected expo- 
nential behaviour of t™ s which, in the asymp- 
totic regime L — > oo, should be given by r SCSD ~ 
exp (2cri 3 ), we perform a x 2 -fit with the ansatz: 



MRS _ rb „cL 



(10) 



which yields Xperdot = 0.897. As a result, we 
find a clear exponential SCSD behaviour for the 
MRS algorihm. On the other hand, for the tun- 
neling times of the MHMC, we expect a monomial 
dependence in L: 



P L q 



(11) 



MRS 



The power law ansatz is well confirmed by the 
fit quality with Xp Cr d o f = 0.795. We also took 
the pessimistic ansatz and tried to detect a po- 
tentially exponential increase of r™ MC . The ex- 
ponential fit gives Xp„ d o f = 0.975. 
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Sustained CPU time in Tflop-hours to generate one 
flip based on tunneling times for MRS (exponential 
fit) and MHMC (lower curve is power law, upper 
curve is exponential fit). The errors of the two ex- 
ponential fits are depicted as dotted lines. The error 
of the power law fit is not visible on this scale. 




In order to compare the efficiency we have to 
take into account the computational effort. The 
complexity of the local Metropolis is given by 
^mrs ~ V whereas the optimized leapfrog scheme 
in MHMC scales as i MH Mc ~ V 5/i [|| . As can be 
seen in the upper figure, the exponential contri- 
bution remains suppressed in the extrapolation. 
A potentially dominating exponential behaviour 
for MHMC can only be detected in future MHMC 
simulations on larger lattices. 

6. Conclusions 

We proposed to make use of the multicanoni- 
cal (MUCA) algorithm within the hybrid Monte 
Carlo (HMC) updating scheme in order to boost 
the tunneling rates. Since both algorithms are 
inherently of global nature, their combination al- 
lows the parallelization of MUCA which could not 
be achieved otherwise. We have demonstrated 
that the fully parallel MHMC algorithm is capa- 
ble to overcome SCSD in compact QED in practi- 
cal simulations, at least up to lattices sizes rs 24 4 . 
On a 24 4 lattice we predict a gain factor of about 
1000 for MHMC over the local metropolis algo- 
rithm with additional reflection. So far, we have 
encouraging experiences on the 18 4 lattice well 
confirming the extrapolations. The investigations 
presented form part of an ongoing study that aims 
at a conclusive FSS analysis of compact QED on 
the Wilson line ||. 
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